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We have achieved the first full implementation of the stochastic projected Gross-Pitaevskii equation for a 
three-dimensional trapped Bose gas at finite temperature. Our work advances previous applications of this the- 
ory, which have only included growth processes, by implementing number-conserving scattering processes. We 
evaluate an analytic expression for the coefficient of the scattering term and compare it to that of the growth 
term in the experimental regime, showing the two coefficients are comparable in size. We give an overview of 
the numerical implementation of the deterministic and stochastic terms for the scattering process, and use simu- 
lations of a condensate excited into a large amplitude breathing mode oscillation to characterize the importance 
of scattering and growth processes in an experimentally accessible regime. We find that in such non-equilibrium 
regimes the scattering can dominate over the growth, leading to qualitatively different system dynamics. In 
particular, the scattering causes the system to rapidly reach thermal equilibrium without greatly depleting the 
condensate, suggesting that it provides a highly coherent energy transfer mechanism. 

PACS numbers: 03.75.Kk, 67.85.De, 05.10.Gg 



I. INTRODUCTION 

Ultracold atomic gases provide a unique environment to 
observe many body quantum phenomena on a mesoscopic 
scale, allowing microscopically derived field theories to be 
readily compared with experiments |[l]|2l. For a Bose gas at 
zero temperature a Bose-Einstein condensate (BEC) forms in 
which essentially all the atoms occupy a single spatial mode 
whose equilibrium and dynamics is well-described by the 
Gross-Pitaevskii equation (GPE) |2, 3|. A growing number 
of experiments have been performed in the finite temperature 
regime, i.e. at temperatures of order the critical temperature 
where many atoms are thermally excited out of the conden- 
sate, e.g. studies of collective modes |i4J, vortex nucleation and 
decay ||5}j7l, condensate growth JS] |9], phase transition dy- 
namics |T|, and superfluid turbulence ifTOl . A quantitative de- 
scription applicable to simulating many of these experiments, 
where thermal fluctuations and dynamics are important, re- 
mains a major technical challenge fTP T^. 

A promising direction of investigation for the finite temper- 
ature regime has been the generalization of Gross-Pitaevskii 
theory to describe the entire low energy part of the sys- 
tem. One approach, typically referred to as the classical field 
method ifTSl IT4ll . involves propagating the GPE with many 
suitably randomized modes. An alternative approach, central 
to the focus of this paper, is to extend the GPE with noise and 
damping terms that represent the coupling to a thermal reser- 
voir of high energy atoms. Recently such stochastic GPEs 
have been applied to a broad range of problems, including de- 
fect formation across phase transitions llTl fTSljTTl . the decay 
of vortices L18r.-20J and solitons 121], and polariton f22| and 
spinor Il23ll24l condensates. The technique has seen quite ex- 
tensive applications to low dimensional systems where ther- 
mal fluctuations can prevent a true condensate from forming 

While phenomenological arguments can be used to obtain a 
generic stochastic GPE, it is possible to derive such a descrip- 
tion from the microscopic theory of a Bose gas. Such for- 
mal derivations have been carried out by the groups of Stoof 



Il231[36l[37] and Gardiner ||38ti40ll validating this approach as 
an ab initio description of non-equilibrium dynamics. 

In this work we present the first complete implementa- 
tion of the stochastic projected Gross-Pitaevskii equation 
(SPGPE), first introduced by Gardiner and Davis in 2003 
||39l . The basic idea of the SPGPE is to sub-divide the sys- 
tem modes into two regions: (i) a high energy region, re- 
ferred to as the incoherent I-region, consisting of sparingly 
occupied modes; (ii) a low energy region, the C-region, of 
appreciably occupied modes, i.e not only the condensate, but 
all other highly Bose-degenerate modes. The use of projec- 
tion operators to define these regions is a key aspect in the 
derivation of the SPGPE theory, and lead to an explicit pro- 
jection operator in the equation of motion for the C-region. 
The thermal effects upon the C-region are described by two 
distinct processes: (i) growth processes where collisions be- 
tween two I-region atoms leads to a change in population of 
the C-region; (ii) scattering processes coiTesponding to colli- 
sion between atoms in the C- and I-regions in which energy 
is transfeiTed but particles are conserved. Simulations includ- 
ing the growth process are relatively straightforward to carry 
out as the noise is additive and uncorrected, and the damping 
term is proportional to the Gross-Pitaevskii evolution opera- 
tor Implementing the scattering process is more technically 
challenging, as the noise is multiplicative and spatially coiTe- 
lated, and the associated damping term involves an intricate 
calculation of cuiTent in the C-region. To date all simulations 
of the SPGPE have only included growth processes |41 J. The 
neglect of the scattering process has been argued as a reason- 
able approximation for systems near equilibrium, but is ex- 
pected to play an important role in non-equilibrium regimes 
02 . Indeed, work on condensate growth within quantum 
kinetic theory JO] |42l |43l has shown that an analogous scat- 
tering reservoir interaction has a large effect on condensation 
dynamics I43ll44l. 

The broad outline of our paper is as follows: 
In Sec. |ll]we briefly review the SPGPE formalism and in- 
troduce the full equation of motion. We also discuss the for- 
mal properties of the various parts of the SPGPE formalism 
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individually, since it has been reasonably common to neglect 
various terms and noises in the equations of motion to arrive 
at simpler theories (e.g. damped GPEs). While the growth 
process gives rise to a grand canonical description of the C- 
region, we show that scattering processes (without growth) re- 
alize a canonical description. By only including the scattering 
process damping term (i.e. neglecting the associated noise) 
we arrive at an energetically damped GPE that evolves any 
initial field to the zero temperature ground state by removing 
energy, but not population. We also evaluate the coefficients 
of the scattering term and compare it to that of the growth 
term in the experimental regime, showing the coefficients to 
be comparable in size. 

In Sec. [Ill] we present our main tests of the formalism, and 
study the evolution of a condensate initial excited into a large 
amplitude breathing mode. We give a number of analytic re- 
sults that we have used to calibrate our numerical algorithm, 
and provide some insight into the effects of scattering pro- 
cesses. Finally we present a study of the breathing mode in an 
experimentally realistic regime and demonstrate that the in- 
clusion of scattering processes causes the system dynamics to 
change in a very significant manner. 

After concluding and surveying the prospects for the full 
SPGPE we provide a detailed overview of our numerical al- 
gorithm in the Appendix. 



II. SPGPE FORMALISM 

Here we briefly overview the formalism of the SPGPE. De- 
tailed derivation of the equation of motion can be found in 
Refs. 



A. The stochastic projected Gross-Pitaevskii equation 

The SPGPE is a c-field method lfT2l where our system is de- 
composed in terms of eigenstates of the single-particle Hamil- 
tonian. 



V{r), 



(1) 



where V{r) — ^m{uj^x'^ + LOyy^ + w^z^), is the harmonic 
confinement potential. The eigenstates of Q, satisfying 
Hsp<l>n{^) — eri0n(r), form a convenient basis of states. Here 
the shorthand n represents all quantum numbers required to 
specify a single-particle state. We introduce a single-particle 
energy cutoff ecut which separates the system into a low en- 
ergy C-region consisting of single-particle modes with eigen- 
values satisfying e„ < ecut- The energy cutoff is chosen so 
that all single-particle modes in the C-region have an appre- 
ciable occupation number, of order unity. In this case we can 
describe the bosonic field for the C-region with a classical 
field lfT2l . ip{r,t), which we represent as a sum over single- 
particle states 



(2) 



where the summation is restricted to all single particle modes 
in the C-region. The remaining high energy incoherent (I) re- 
gion contains thermally occupied modes, and acts as a thermal 
reservoir for the C-region. 

The I-region is assumed to be in equilibrium with a well de- 
fined temperature T and chemical potential /i, and in a semi- 
classical local-density treatment is described by the single par- 
ticle Wigner function 



F(r.k) = 



1 



exp[(Sw(r,k)-Ai)/fci3r]-l' 



(3) 



where the energy is huj{r, k) = S^k^/2m + V{r). 

Accounting for the reservoir interactions with the C-region 
leads to the non-local Stratonovich stochastic equation of mo- 
tion for tp{r, t), known as the SPGPE (writing ip = ^(r, t) for 
brevity) | 



(S) = dij\^ + dVlc + {S) dt/j] 



M' 



where 
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in- C)^pdt + dWG{r,t) 



VM{r, t)'ipdt + ii^dWuir, t 



(4) 

(5) 
(6) 
(7) 



with (5) denoting Stratonovich integration 11451 . In 
Sees. |II A2prA3| we define the terms given in Eqs. (|5]l-(|7]i, 
respectively. A common feature to all terms is the projection 
operator defined by 

P/(r)^^<^„(r) [ dVC(r')/(r'), (8) 
which formally restricts the evolution of to the C-region. 



1. Hamihonian evolution term: dV'j^ 

The Hamiltonian evolution operator for the C-region, £, is 
the usual GPE operator 



(9) 



where u — ATrfi^a/m with a the s-wave scattering length. 
However, as it appears projected in Eq. (|5]l this equation by 
itself is refen-ed to as the projected GPE (PGPE). The PGPE 
contains the important interactions between the low energy 
modes as schematically indicated in Fig.[TJa). 



2. Growth term: dip | ^ 

The growth process [Eq. (|6|] describes the collisional in- 
teraction in which two I-region atoms collide leading to pop- 
ulation growth of the C-region [see Fig.[TJb)]. It is set by the 
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FIG. 1. The different interparticle scattering processes arising from 
each term of the SPGPE (j4]l, where ecut specifies the boundary be- 
tween C and I. The I-region is static, parameterized by the tempera- 
ture (T) and chemical potential /i. (a) Non-linear mixing interaction 
between modes in the C-region. (b) Growth process in which both 
energy and number are transferred between C and I. (c) Scattering 
process in which energy is transferred between C and I. 



the functions G(r) and M{r) (providing G{r) ^ 0), the 
SPGPE evolves arbitrary initial conditions towards samples 
of the grand canonical ensemble, with equilibrium probability 
P(V') oc cxp{-Kc/kBT), where Kq ^ Eq - fiNc is the 
grand canonical Hamiltonian, with 

Ec = J rHsp^ + 1 / i^^i^' ^^^'^ 

Nc^ I d^r IVP, (15) 



the C-region energy and atom number, respectively. 

In the next three subsections we examine the various sub- 
theories of the full SPGPE obtained by neglected various 
terms. This gives some insight into the role of the individual 
terms. Some of these sub-theories have been quite extensively 



used in the field, particularly the established PGPE [Sec. |IIC | 
and simple growth SPGPE [Sec 



growth rate G(r) (see Sec. InDl) and the Gaussian complex sub-theory presented in Sec.|irE 



TTdI . The scattering SPGPE 



noise dWciic, t) has the non-zero correlation 

{dW^ir, t)dWG{r', t)) = 2G(r)Jc(r, r')dt, (10) 

where (5c(r,r') = J2nec '^^i-(^)'f'n(^') ^ delta function in 
the C-region. 

3. Scattering term: dip | 

The scattering process [Eq. (j?])] transfers energy and mo- 
mentum between the C and I-region without population trans- 
fer [see Fig.[T|c)]. It is set by the rate Af(r) (see Sec. IIEi 
which determines the effective potential 

VM{r) = J dh' M(r - r')V' • j(r'), (11) 

that couples the divergence of the C-field current 

ih 



^has not been considered be- 
fore. In Table|l]we present a summary overview of the various 
sub-theories considered and their general properties. 



C. Projected GPE 

The PGPE, defined by Eq. (|5]), is formally (and numeri- 
cally) number and energy conserving, for any finite cutoff ecut 
and single-particle basis. Since VV = V, we have Vip = ip 
and 



dNr 



dt 



H 



d^r tP*V 
d^r {V'4j) 



':£tpdt 



-Ctpdt 
n 



h.c. 
+ h.c. 



= 



(16) 

(17) 
(18) 



j(r) 



2m 



(12) 



to the I-region. The scattering noise dWA/(r,t) is real, de- 
fined by its non-vanishing correlation 



{dWM{v,t)dWM{v',t)) = 2M(r - T')dt. 



(13) 



The scattering noise forms a real-valued stochastic potential 
and thus generates a phase-diffusion process for the C-field 
evolution that transfers energy and momentum between the C 
and I-region atoms while preserving their populations. 



B. Formal properties of the full SPGPE 

Here we briefly overview the formal properties of the full 
SPGPE theory. Due to the coupling of the C-region with 
the reservoir, the SPGPE [Eq. gives a grand canoni- 

cal description of the C-region. Irrespective of the form of 



where to get (17 1 we have used the fact that V is Hermi- 
tian ll46l . Similar (lengthier) reasoning gives dEc/dt\^ = 0. 
While a formally Hamiltonian theory, the nonlinear interac- 
tions generate ergodic dynamics and the equation samples the 
microcanonical ensemble in equilibrium [47 1. The PGPE de- 
scribes the dynamics of both thermal and coherent C-region 
atoms non-perterbatively, giving a quantitatively accurate de- 
scription of finite temperature systems in (or near) equilibrium 
where the reservoir interaction can be neglected |48|. Hence 
the PGPE has mainly been applied to the equilibrium proper- 
ties of the finite temperature Bose gas [13. .48-58] . Dynamical 
studies have considered vortex nucleation ||59l and collective 
modes [60l. 

Numerically, the PGPE is a fully dealiased spectral method 
for propagating the GPE. In the basis of plane waves (i.e. the 
Fourier spectral method), the PGPE eliminates all spurious 
aliasing generated by four-wave interactions. In practice this 
is achieved by evaluating the interaction term for a wave func- 
tion of 71 points per spatial dimension {n momentum modes) 
on a grid with 2n points in each dimension (extending out to 
twice the momentum cutoff of the wavefuction) [49|. This 
procedure is easily generalized to other bases IJ5u50J . 
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Method 


Ensemble 


Conserved quantities 


Variable quantities 


PGPE[Eq. (5i] 


Microcanonical 


Ec,Nc 


- 


Simple growth SPGPE [Eq. 


20 ] 


quiet 


NA 


- 


Ec,Nc 


noisy 


Grand-canonical 




Ec,Nc 


Scattering SPGPE [Eq. |27f] 


quiet 


NA 


Nc 


Ec 


noisy 


Canonical 


Nc 


Ec 


(full) SPGPE [Eq. ( 


4|, 


Grand-canonical 




Ec,Nc 



TABLE I. Summary of the different theories considered in this paper. Quiet implies that the appropriate noise term is neglected, giving 
the damped PGPE (quiet simple growth SPGPE) and energetically damped PGPE (quiet scattering SPGPE). The different reservoir processes 
leading to the distinct ensemble descriptions are summarized in Fig.^ We emphasize that the scattering SPGPE enables a dynamical canonical 
description of the c-field. The c-field methods which are implemented for the first time in this work are shown in bold red font. 



D. (Simple) Growth SPGPE 



3. Properties of the simple growth SPGPE 



The only form of the SPGPE Q which has been used for 
numerical simulations is the simple growth SPGPE, in which 
scattering processes are neglected and the growth rate G(r) 
is taken as spatially uniform. The resulting equation is easily 
handled numerically and closely connected to the Ginzburg- 
Landau (j)'^ theory. 



1. Growth rate 

When the I-region is near equilibrium and described us- 
ing Eq. ([3]l, G(r) is approximately spatially constant over the 
condensate 1 15 1. In this case the growth amplitude can be cal- 
culated explicitly as 1151 



G(r) « 7 



70 5] 



(19) 



where 70 = Sa^/A^^, with A^s = yj2TTh? /mkBT, /3 = 
l/ksT and u, w] the Lerch transcendent. 



Without noise (dW-^ = 0, a case we call quiet), the sim- 
ple growth SPGPE reduces to the damped PGPE 1631 . The 
damped PGPE evolves Nq to equilibrium: 



dNc 
dt 



H-\-'y, quiet 



- -27[/l(i) -ii]Nc, 



(23) 



where fi{t) = J d'^rip* £ip/Nc is the instantaneous chemi- 
cal potential. This evolution also causes energy to decay uni- 
formly: 



dKr 



dt 



H~\-'y , quiet 



2^7 



d^v\{^Ji~C)4'\\ (24) 



thus minimizing Kq and damping out thermal fluctuations. 
The equilibrium solution is the zero temperature ground state 
of the PGPE ^ satisfying ^iipa = 7'{/:i/'o}- 

Upon reintroducing the noise, Eq. pO] ) samples the grand 
canonical ensemble and ifciV"] > ^c[^o] for any sample 
V^. However, all equilibrium properties are independent of the 
choice of 7. Confirming this property provides an excellent 
test of any numerical implementation. 



Evolution equation 



E. Scattering SPGPE 



The simple growth SPGPE is given by 



dm 



H+1 



dm 



H 



dm 



where 



dtp\ 



1 



(fi- C)'ijjdt + dW-^{r,t) 



and the noise correlation is 

{dW;ir,t)dW^{r',t)) = 2-fScir,r')dt. 



(20) 



(21) 



(22) 



The numerical integration of (j2T]i is of a similar computa- 
tional expense to solving the PGPE |61 1 since the noise is ad- 
ditive and weak, and we can use a higher order Runge-Kutta 
algorithm to achieve stochastic convergence 1621 . 



We define the scattering SPGPE to be the sub-theory of the 
SPGPE obtained by neglecting growth terms. 



1. Scattering rate 
The scattering rate is given by 

M(r 

where 



(2^) 



zk-r 



M 



m-Ka^keT 1 



(25) 



(26) 
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In Ref. 1391 this rate was referred to as the "simplified non- 
local form", however it can be shown (see Appendix [A]) that, 
within the semiclassical approximation used in the Wigner 
function F(r, k) description of the I-region, it is exactly equal 
to the full non-local form. 



Evolution equation 



The scattering SPGPE takes the form 



(5) d^p\ 



H+M 



H 



(S) d^\ 



Af ■ 



(27) 



As the scattering process is fundamentally number- 
conserving, the scattering SPGPE provides a dynamical 
canonical description of the C-region. Equilibrium states 
of ( [27] i sample the canonical distribution with probability 
P{^P) cx expi^Ec/kBT). 

Numerically integrating the scattering SPGPE ( p7| poses a 
heavier technical challenge compared to evaluating the sim- 
ple growth SPGPE (20 1. Due to the multiplicative noise aris- 
ing from the scattering process, we are restricted to using a 
semi-implicit Euler algorithm to evolve the stochastic differ- 
ential equation (see appendix section B 2 1. This algorithm 
is first order in the weak sense of convergence of Ref. B31 . 
hence is much more inefficient (in computation resources) 
than the higher order Runge-Kutta algorithms that can be used 
to evolve the simple growth SPGPE. Further complexity arises 
through the non-local deterministic term, VmIi") [Eq. (Ill], 
and in sampling the correlated noise dWM- The key steps in 
our numerical implementation use techniques developed for 
solving the PGPE with dipole-dipole interactions |64], and a 
full overview of the algorithm is given in Appendix [B] 



3. Properties of the scattering SPGPE 



Setting dWu = in Eq. ( 27 1 leads to the quiet form of 
the scattering SPGPE, which we call the energetically damped 
PGPE. The deterministic part of the scattering term Vm takes 
the form of an effective potential that acts to reduce the energy 
of the C-region, Eq. ( 14i. Using Eq. (27 1 with dWM = we 
find 

dKc 



dt 



dEc 



H+M, quiet 



dt 



H+M,quiet 



J d3rFM(r)V-j(r) 



Substituting Eq. (Ill and Eq. ( 25 1 gives 



dKc 



dt 



H+M,quiet 



(28) 



(29) 



which is negative semi-definite. The scattering term causes 
the C-region energy to monotonically decrease, which pro- 
vides a useful consistency check on the the accuracy of the 
numerical evaluation of Vm- 



The Ehrenfest relation ( 28 i gives an important physical in- 
sight into the role of the scattering process in the SPGPE. 
The scattering process generates a dissipative interaction that 



enters the evolution as an effective Hamiltonian term, in the 
form of a stochastic effective potential. When the noise is 
neglected, energy is continually removed through evolution 



[Eq. ( 29 1] and the system proceeds toward the absolute ground 



state. The inclusion of the noise ensures that the effective po- 
tential is stochastic, maintaining the finite temperature char- 
acter of the C-region. 

We have identified some basic tests that proved useful in 
validating our numerics. First, equilibrium ensemble prop- 
erties will be independent of the scattering coefficient M. 
for a given temperature (much as the simple growth SPGPE 
equilibria are independent of 7), and should be equivalent to 
those generated by the simple growth SPGPE for the same fi- 
nal particle number (assuming equivalence of the canonical 
and grand canonical ensembles). Second, we expect that a 
stochastic scattering term implementation will evolve Nq par- 
ticles in the C-region into thermal equilibrium. Evolution ac- 
cording to the deterministic term will evolve the c-field region 
toward the A^c-particle PGPE ground state. 



F. Comparison of the scattering and growtli coefficients 

To assess the relative importance of scattering and growth 
terms in the SPGPE it is of interest to evaluate the respective 
coefficients (70 and M.) in regimes relevant to experiments. It 
is useful to consider these in their dimensionless forms 



7o 



M 



Mh 

kuTxl 



xt IJ- 



(30) 
(31) 



where xq = \JhJmLj with tj = ^fnT^^J^^. The approxima- 
tion in Eq. ( [3l] l is obtained for the usual validity regime of 
the SPGPE theory [12| kuT > ^, and typical cutoff choice 
Ecut ^ 3/^. Thus the ratio of the coefficients is given by 



M 

70 



(32) 



In usual experimental regimes A^s ^ (for temperatures 
near the critical temperature) and thus we conclude that the 
scattering coefficient is significant, potentially appreciably 
exceeding the growth coefficient. 

To be more quantitative we evaluate the coefficients using a 
numerical calculation (see Ref. 1 18|) for a spherically trapped 
system with 0;^ = 27r x 10 Hz and a range of total atom num- 
bers N = iVc + A'l (where N\ is the number of atoms in the 
I-region) and temperatures (^ is found to ensure N is fixed as 
T varies). The results for 7 and M. are shown in Fig. [2] noting 
that for the growth term we evaluate the full coefficient 7 ( [T9| . 
These results support the qualitative analysis given above, and 
show the coefficients are similar in size over a broad regime. 
However the net effect of scattering depends on the diver- 
gence of currents in the C-region, and can be small for quasi- 
equilibrium dynamics where the simple growth SPGPE has 
been successful 1201 . 
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X 10 



N = 1 X 10* 
K = 1 X 10"' 
Af = 1 X 10° 




FIG. 2. The growth (7) and scattering (A4) coefficients, and their 
ratio, as a function of /i/kBT. Reservoir parameters T,fi, tcut are 
determined from the Hatree-Fock parameter estimation scheme from 
Ref. |18|. The parameters found ensure T can be varied without 
altering the total atom number A'^, and critical temperature Tc{N). 
Along each curve the temperature varies from T — TctoT — OATc, 
where Tc corresponds to a total atom number of A'^ — 1x10* (black), 
TV = 1 X 10^ (dashed gray), iV = 1 x lO'^ (gray). 



III. BREATHING MODE DECAY 

As an application of our implementation of the SPGPE we 
study the case of a BEC excited into a large amplitude breath- 
ing oscillation. 



field, found analytically by evaluating ( 14 1, is 



Nc, 



4to 



1 



2 2 
an 



3ma;"(T 



2^2 



4 V327r3/2CT3 

(34) 

We begin by looking at a spherically trapped condensate, 
with a trapping frequency of = 27r x 10 Hz, and choose 
parameters of the initial field to be = xq, k 

Nc,^ = 1 



Xq , with 

X 10"' ^^Rb atoms. For these parameters, we find 
from (|34|l that Eq = Ec{t = 0) = 8.44to7Vc,i- The ground 



state of the Gross-Pitaevskii equation for this system is Eq = 
3.53hujNc.i, so our initial state is far from equilibrium. The 
current density of this initial condition is given by 



j(r) = tkIV'I' 



(35) 



i.e. a radially expanding motion. Using (35 1 we can evaluate 
( 1 1 1 to find an analytic expression for Vm 



g{r/a)huj, 



(36) 



where 



V2a 



V2x] e-''\TfL{x), (37) 



with erfi(a;) = — ierf(ia;). For x <C 1, g{x) has the asymp- 
totic expansion 



- 11- -X' 



+ 0(x4), 



(38) 



while lim^_j.oo g{x) — 0. For our choice of initial condi- 
tion (k > 0, i.e. radial flow away from the origin) we have 
Vm{0) < and VA/(r) — > as r — >■ cx). Thus for small r 
Vm (r) is approximately described as an additional harmonic 
potential well which acts against the radial expansion. In ad- 
dition to providing extra confinement, V^/ has a negative en- 
ergy offset, which modifies the effective energy minimum of 
the potential experienced by the system, and thus could effect 
growth into the condensate. 

In Fig. [3] (upper right subfigure) we compare Vj(\{r) with 
VM{r) obtained numerically using the initial condition (33 1 
and find excellent agreement. 



B. Scattering: deterministic versus noisy dynamics 



A. Gaussian wave function 

For the results we consider in this section we begin with the 
well-defined initial condition for the C-region field: 



(33) 



where iVc ^ is the initial C-region number, and a and k are 
real constants. The C-region energy per particle of this initial 



We examine the scattering SPGPE evolution of a Gaussian 



wave function (33 1, with the parameters specified in section 



III A First, we consider the effect of the scattering effective 



potential term on its own by using the energetically damped 
PGPE [Eq. (27 1 without noise]. In Fig. [3] we show the deter- 
ministic evolution of the Gaussian wave function, comparing 
the density with the effective potential at a range of repre- 
sentative times. The arrows on the density images show the 
direction of the breathing motion (i.e. the direction of current 
flow) at each time. As the breathing motion evolves, VmI?") 
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Column Density 



Potential 




-4 -2 2 4 

X (xo) 



FIG. 3. Results of the energetically damped PGPE, comparing the 
density with the scattering effective potential where M — 0.005, at 
representative times (in units of trap cycles). Column densities are 
shown in the left column, with arrows indicating the direction of the 
breathing motion and flow of current at each time. The scattering po- 
tential, Vm (red curve), and a radial slice of density in energy units, 
the local C-region interaction energy it|i/'P (black curve), are shown 
in the right column. At t = we show Vj(j (blue curve), the ana- 
lytical scattering potential for a Gaussian wavefunction found from 
1 36 1. (blue dashed curve), is i36i found using the harmonic 

approximation in Eq. OSl. 



acts against the density change. As the condensate expands 
we have VA/(r 0) < 0, where the negative potential acts 
against radial expansion. Then as the condensate begins to 
contract and the flow is directed towards the origin, we see 
VM{r ~ 0) > (see t = 0.39 eye). This is consistent with 
Eq. p6l ) evaluated with k < (i.e. consistent with an inward 
current). Finally the currents are completely damped out and 
the system reaches equilibrium, i.e. the ground state of the 
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FIG. 4. C-field energy per particle as a function of time (in trap cy- 
cles), for the scattering SPGPE evolution of an initial Gaussian wave 
function with A'c = 1 x 10** atoms, (solid lines) scattering SPGPE 
evolution including the noise, at temperature of T = IGhuj/ks for a 
range of M . (dashed lines) Evolution with the deterministic scatter- 
ing term, neglecting the noise, for a range of M . The initial energy 
(-Eq) determined analytically from \3A) agrees with the numerical 
calculation. For the quiet simulations, the ground state energy agrees 
with the Gross-Pitaevskii equation ground state energy labelled Eg- 
For the noisy simulations, the occupation of thermally excited C- 
region modes raises the equilibrium energy. 



PGPE (shown at i = 5 trap cycles). 

The evolution of the C-region energy for the energetically 
damped PGPE is shown in Fig.|4] We see the expected behav- 
ior, namely that Vm acts to reduce the energy monotonically 
until the system reaches the ground state, with energy con- 
sistent with the zero temperature GPE ground state {Eq). The 
final state is independent of the value of M., although this does 
influence the rate at which the final states are reached. We can 
quantify the initial effect that Vm('') has on the rate of change 



of the C-region energy. To do this we evaluate Eq. ( 29 1 for the 
initial field ([33]l, giving 



dEc 
dt 



(39) 



which describes the rapid loss of energy initially seen in 
Fig. |4] Using the parameters of our initial state, we find 
dEc/dt — —5.1 X IO^/kD^, which agrees with our numeri- 



cal evaluation of Eq. \29\ to better than one part in 10"*. 

We now examine the role of the scattering noise term, for 
the same parameters as above, and a temperature of T = 
lOhui/kB- The energetic evolution for this case is also shown 
in Fig. |4] For the stochastic simulations, we have averaged 
over 10 trajectories for each parameter set. From the initial 
condition we again observe energy to decay, however unlike 
the noiseless simulations where it strictly decreases, we see 
that this is not the case for the noisy simulation. The local 
peaks in Eq occur when the condensate is fully contracted 
in its breathing motion. We observe that, as expected, the fi- 
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nal finite temperature equilibrium state is independent of A4. 
The equilibrium states has appreciably more energy than the 
ground state, reflecting the thermal excitation of the system. 



X 10 



C. Finite temperature breatliing mode decay in 
experimentally realistic regime 



an 



Finally we extend our study of the breathing mode to a 
regime with experimentally realizable parameters and com- 
pare the predictions of the full, scattering, and simple growth 
SPGPEs. 



1. Parameter choice 

In order to give a well defined comparison, we choose phys- 
ically consistent reservoir parameters using the Hartree-Fock 
parameter estimation scheme described in Ref. ifTSl . For a to- 
tal atom number of iV = A'c.i + A'l = 5 x 10* atoms, we find 
T = 29.4to/fcB, = 4.8to, Ccut = 15.9to (T sa O-SSTJ, 
giving (7, M) = (1.5, 2.7) x 10"'' [from Eqs. ^ and E^]. 



Our initial state consists of a Gaussian field with radial phase 
gradi ent ( [33] l, with the same parameters as specified in section 
III AI except with Na — 1.22 x 10*, which corresponds to 



the Thomas-Fermi condensate number from our value of /i. 
We evolve this initial state with the same reservoir parameters 
T, jjL, and A'l determined above. 



2. Comparison of evolution 

In Fig. |5] we compare the evolution of A^c and Kq as an 
average of 50 trajectories, throughout the decay of the breath- 
ing motion. The number and energy reach equilibrium much 
faster for the full and scattering SPGPEs than for the simple 
growth SPGPE. The behavior of Kq is almost identical for the 
full and scattering SPGPEs, with the difference in the value of 
Nq that these two theories equilibrate to arising because the 
scattering SPGPE conserves number In contrast, the simple 
growth SPGPE predicts very different behavior: Kq changes 
in a similar way to the other theories, but on a much slower 
timescale. Nf^ instead evolves in a different manner, decreas- 
ing to about 80% of Nc.i before slowly returning towards this 
initial value (an equilibrium value of Nq, = 1.15 x 10* is 
eventually reached). Note we have verified that our numeri- 
cal algorithm produces the same equilibrium state (after suf- 
ficiently long times) for the full and simple growth SPGPEs, 
irrespective of the value of 7 and M.. 

To quantify the decay of the actual breathing mode we cal- 
culate {r"^), which provides a measure of the average system 
width (noting (r) w 0). The evolution of (r^) is shown in 
Fig. [5] with the initial oscillations of (r^) correspond to the 
breathing mode of the condensate. We see the damping of (r^) 
for all methods is broadly consistent with that of Kq. How- 
ever (r^) for the simple growth SPGPE simulations shows an 
interesting difference: The oscillations decay within 10 trap 
cycles, despite (r^) (as well as Ac and Kq) being far from 
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FIG. 5. A^c, Kc, and (r^) as a function of time, for the evolution of a 
Gaussian wave function using the SPGPE (black), scattering SPGPE 
(red), and simple growth SPGPE (blue). We show the time evolution 
for each method until equilibrium has been reached. In the inset 
we focus on the early time dynamics, over which the SPGPE and 
scattering SPGPE reach equilibrium at a much faster rate than the 
simple growth SPGPE. The SPGPE and scattering SPGPE results 
for Kc and (r^) are virtually indistinguishable. 



equilibrium. After this time (r^) decays slowly (without os- 
cillation) towards the equilibrium value. 



3. Condensate dynamics 

To understand the marked difference between the simple 
growth evolution and the other two theories it is useful to con- 
sider the behavior of the condensate, which we examine in 
Fig. |6] We determine the condensate number from C-region 
field using the Penrose-Onsager definition [651: We form the 
one-body density matrix 



Pi(r,r',t) = (V'*(r,i)^A(r',t)), 



(40) 
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FIG. 6. Condensate fraction found from the Penrose-Onsager cri- 
terion, as a function of time for tlie same system as in Fig. |5] The 
SPGPE (blaclc) and scattering SPGPE (red) simulations are shown in 
the inset, and the simple growth (blue) simulations are shown in the 
main figure. 



where angle-brackets denote an ensemble average over tra- 
jectories at a given time t. The condensate number No{t) is 
defined as the largest eigenvalue of pi (r, r', t), and in our cal- 
culations always greatly exceeds the next largest eigenvalue. 
Our initial state for the C-region is essentially a pure conden- 
sate, with A'o = A^c at t = 0. Similar to the observations 
made of Fig. |5]we see that Nq reaches equilibrium far more 
rapidly for the full and scattering SPGPEs compared with the 
simple growth SPGPE. Our results show that in the simple 
growth SPGPE the condensate fraction rapidly drops (over a 
time period consistent with the rapid decay in (r^)) to a min- 
imum condensate number of Nq — 1.85 x 10'^ at t — 8 cy- 
cles. The very slow approach to equilibrium observed after 
this time corresponds to a re-condensation process, as can be 
seen in the long-time evolution of (r^) in Fig.jsj^c). The scat- 
tering term (i.e. full and scattering SPGPEs) allows a different 
and very effective route to rapidly dissipate the energy of the 
breathing mode without drastically reducing the condensate 
fraction. 



IV. CONCLUSIONS 

When the full SPGPE theory was derived by Gardiner and 
Davis in 2003 they touted that "This approach is distinguished 
by the control of the approximations made in its derivation, 
and by the feasibility of its numerical implementation" ll39l . 
In this work, we have finally realized a numerical implemen- 
tation of this theory, and demonstrated practical simulations 
in the experimental regime. 

To date all applications of the SPGPE have been made 
within the simple growth approximation in which the scat- 
tering terms are neglected. Using our algorithm we are able 
to assess the effects of the scattering terms. We have verified 



that, when growth terms are neglected, the scattering terms 
evolve the system to an equilibrium state that is independent 
of the scattering amplitude coefficient (M), and that samples 
the canonical ensemble for the C-field region. The latter prop- 
erty is distinct to the simple growth and full SPGPE descrip- 
tions that exchange both energy and particles with the reser- 
voir, and hence sample the grand canonical ensemble in equi- 
librium. 

We have applied our theory to study the evolution of a finite 
temperature condensate excited into a large amplitude breath- 
ing mode in a physically realizable regime. An important, and 
somewhat surprising observation, is that the SPGPE with the 
scattering terms predicts a qualitatively different evolution to 
the simple growth SPGPE: with the inclusion of scattering, the 
breathing oscillation is efficiently damped without greatly de- 
pleting the condensate, allowing equilibrium to be established 
on a much shorter time scale. The energy damping is due 
to the scattering effective potential that precisely opposes su- 
perfluid motion, and the results suggest that the scattering de- 
scribes coherent energy exchange with the reservoir, a striking 
consequence of Bose-enhancement. Our results indicate that 
scattering terms may be important in highly non-equilibrium 
regimes encountered in ultra-cold gases, and that quantitative 
evidence for the dominance of scattering over growth might 
be easily measured in experiments. 

Future work with the full SPGPE will be to advance our 
understanding of non-equilibrium dynamics in the finite tem- 
perature regime. An exciting prospect is the direct compar- 
ison with experiments of a non-equilibrium scenario of fi- 
nite temperature dynamics, such as the breathing mode decay 
studied here, or the dynamics of condensate growth during a 
quench 
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Appendix A: Scattering Rate 

The full non-local scattering rate, as derived in lf39ll is of 
the form 



M(R,r) = 



1 



(Al) 



where 



£'nun(R) = max Ecut , 



8m 



14ff(R) . (A2) 
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This result is obtained as the Fourier transform of 



M(R,k) 



47ru2 rd^Ki fd^K 
<5K-c^2)m,Ki)[l + F(R,K2)], (A3) 



^(5(Ki -K2-k) 



where Ki and K2 are the wavevectors of the I-region 
atom before and after a scattering collision, respectively [see 
Fig. [TJc)]. The I subscripts on the integrals in Eq. (A3i in- 
dicate the domain of integration is restricted to the I-region, 
where /ki;(R, k) > e^ut- 

This rate can be further simplified as follows: In a semiclas- 
sical description conservation of momentum during the scat- 
tering collision requires Ki + K3 = K2 + K4, where K3 
and K4 represent the before and after wave-vectors for the C- 
region atoms participating in the collision. By definition the 
C-region atoms must satisfy 



2m 



Kff(R)<ecut, J = 3, 4. 



(A4) 



Since the momentum transfer in the scattering event satisfies 
?ik = KKi — fiK2 = /1K4 — /1K3, we have 



;i2|k|- 

8m 



1 



K=ff(R) = ^ 



2m 

+Kff(R) 



1 

< 

- 4 



2m 

+Kff(R) 

< Ecut- 



+ 



+ 



2m 



2m 



(A5) 



/i2|K,||K4 



(A6) 
(A7) 



Thus £',„in(R) = Ecut (in Ref [39| this was stated as an 
approximation), and the R dependence is lost in Eq. ( |Al| l, 
i.e. M(R, r) — > M (r), which is the form we use in this work. 



Appendix B: Outline of the numerical algorithm for the 
scattering SPGPE 

Here we detail how we efficiently implement the scatter- 
ing SPGPE in the harmonic oscillator basis. The numerical 
implementation of the simple growth SPGPE has been previ- 
ously outlined |fl2l|T5|, being only slightly more complicated 
than integrating the PGPE for a harmonically trapped system 
W^- Thus we focus here on our evaluation of the determinis- 
tic scattering effective potential and the associated scattering 
noise. We emphasize that here we present a simple overview 
of the method and a full and detailed account of our algo- 
rithm, particularly the use of quadratures to accurately and ef- 
ficiently evaluate the necessary matrix elements, wiU be given 
elsewhere. 



in terms of the basis of harmonic oscillator modes </)„ of the 
single-particle Hamiltonian satisfying -ffsp^n — ^n4>n, where 
c„ are time dependent complex amplitudes and n represents 
all quantum numbers required to specify a single-particle 
state. This choice is convenient because it allows us to ef- 
ficiently implement the projection by restricting the spectral 
modes [as indicated in Eq. ( |Bl| l] to the set 

C = {n : e„ < ecut} • (B2) 



Projecting the scattering SPGPE ( [27] i onto the basis-set 
modes in the C-region we obtain a system of equations for 
the evolution of the amplitudes, i.e. 

(5) dcn = -i[e,iC„ + Gn + Sn]dt + ^ Bn,n dw,n, (B3) 



where 



G„ 



u / d3r0;(r)|V(r,i)|V(r,i), 



S„ = I d\<t>l{v)VM{v,t)i^{v,t), 
I / d^r4>l{r)iP{r,t)xmir) 



(B4) 



(B5) 



(B6) 



where we introduce the functions Xm(r) later [see Eq. (B25 1] 
and dwm is the standard real Wiener process satisfying 
(dw„) = 0, {dw.^dwn) = Smndt. 

There are two main steps in solving this equation: (i) time- 
evolution to step this equation forward in time; and (ii) eval- 



uating the non-linear matrix elements ( B4 1-( B6 1 at each time 
step. 

2. Time-evolution 

We employ the weak vector semi-implicit Euler algorithm 
||45, 62 67, 68] to evolve our stochastic equations forward in 
time. As this algorithm is extensively discussed in the litera- 



ture we just briefly review the algorithm here. Equation (B3 1 
is of the general form 

(5) dCn = a„(t, c)dt + ^ S„„(i, c)(iw„(i), (B7) 
j 

where a„ = — i [e„ c„ + G„ + S'n ] , and we use the notation c to 
represent the dependence of matrix elements on the full field 
(V'). The solution is propagated to a set of discrete times tj = 
j At, where At is the step size, and we denote that solution 
at time tj as Cn \ Using this solution, the solution at the next 
time-step is computed as Cn^^'' ~ cP^ + Ac'^^\ where 

Ac(f' ) = a„ (t, , c(^") ) Ai + ^ (tj , c(^) ) Aw^J^'> , (B8) 



with 



1. Basis state representation 

We use a spectral representation of the c-field 

nec 



(Bl) 



— 2 ' n J-> 

- _ 1 



(B9) 

(BIO) 
(BID 
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Note formally Awn ^ = J^'*^ dwn, however in practice we 



sample Aw^J as a real Gaussian distributed random variable 
of variance At [c.f. Eq JbTTJi]. 



Because Cn in Eq. (B8 1 depends on c„ ' [i.e. (|B9|)], this 



equation is implicit. In practice a few iterations of Eq. ( B8 1 are 
usually sufficient to obtain convergence at small step-sizes. 

This algorithm is correct to 0{At^) in the limit of zero 
stochastic noise, and is convergent for the stochastic problem 
with a strong order of At^/^ (i.e. for individual trajectories) 
and with a weak order of At^ (i.e. for quantities calculated in 
the distribution). 



c. Scattering noise 



While the scattering noise ( 13 1 has a non-local correlation 
function in position space, in Fourier space it satisfies 

{dWM{k)dWM{k')) = / d^r / d'r' g-"''--' 



(27r)3 , 
X / d^q 



2Mdt 



„iq-(r-r') 

(5(k + k'). 



(B15) 
(B16) 



3. Matrix elements 

As noted above the matrix elements of the usual Gross- 
Pitaevskii evolution [i.e. dB4|] and stochastic growth are dealt 
with elsewhere lfT2l [TSl 16 111 and we do not repeat this here. 
Instead we focus on the two new terms associated with scat- 
tering. 

a. Use of Fourier transforms to simplify scattering terms 
In what follows we will use the notation 



d^r 

(27r)3/2' 



7(r), 



(B12) 



to denote the three-dimensional Fourier transform of the func- 
tion /(r), which could be either a scalar or vector function, 
with associated inverse transform /(r) — T~^{f(k)}. 

These transforms can be efficiently and accurately imple- 
mented in the basis-set approach using the fact that the oscil- 
lator basis is the eigenbasis of the Fourier transform operator 
We note that the Fourier transform of the field (e.g. J^{i/;}) and 
density-type quantities (e.g. J^{|?/;p}) are treated in different 
ways, as discussed in Ref. L64J . 



b. Scattering effective potential 



That is, the noise is anti-diagonal in k-space. We implement 
the noise generation procedure using the following result: 

Choosing the oscillator basis to be separated into prod- 
ucts of one-dimensional oscillator eigenstates (/i,i(r) — 

Vri} {x)tpny {y)^n} (z), whcrc wc havc decomposed the quan- 
tum number as n = {rix, ny,nz}- These ID states are even or 
odd, as determined by their quantum number, e.g. tp^i] (x) — 
(^_X)"^iy9l^-'(— a;), and thus full basis states have the parity 
property 



Mr) - i-irM-r), 



(B17) 



where cr — + riy + Uz- Taking the oscillator basis to be 
real (in position space), the Fourier transformed basis modes 
(/)„(k) are then purely real or imaginary functions depending 
on this symmetry, i.e. 



i„(k) ^ (-z)^$„(k), 



(B18) 



where the $,i(k) are a purely-real set of orthonormal or- 
bitals. Importantly, the functions <l>„(k) are, to within a scal- 
ing along each dimension, identical to the position space func- 
tions (j)n{r), and thus property (B17 1 holds 



$„(k) = (-l)^$„(-k) 



(B19) 



Using this result, the scattering noise is then constructed in 
k-space as 



Using ( 25 I and the convolution theorem, the scattering ef- 
fective potential can be evaluated as 



VM{r) 



^-i{ik. J-{j(r)}}, (B13) 



where k — k/|k|. The current j(r) is quite conveniently 
formed in the oscillator basis making use of the step opera- 
tors to take the spatial derivatives. 



Equation ( B 1 3 1 reveals that the scattering effective poten- 
tial depends on the radial part of the current in k-space. This 
corresponds to the irrotational part of the of the current, i.e. 



(B14) 



where j(r) = j||(r) + j_L(r) is the Helmholtz decomposition 
withVxj||(r) = V-jj_(r)=0. 



rfW^(k) = J^^^„(k)dw;„, (B20) 



which has the desired correlation function ( |B16| l: 

2Mdt 

2Mdt 



dW{k)dW{k.')) = ^^^(-l)'^$„(k)$„(k'), 

n 

(B21) 

^$„(k)$„(-k'), (B22) 



|k| 
2Mdt 



'5c(k,-k'), 



(B23) 



where 5c(k, k') = E„ <fn(k)$„(k') = E„ C(k)<^„(k') is 
the projected delta function. 
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Thus, in r-space the noise is given by 

dWM{r, t)=Y^ Xm(r) dw^, (B24) 



where we have introduced 



X„^(r)=^|./^0„(k)l. (B25) 
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